﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace PhaseUatt3
{
    public class regulatorConditions
    {
        public double[] w;
        public double[] f;
        public double[] df;

        public int maxStep;
        public double[] h;

        public double[] s; 

        public double m; //масса шаров центробежного регулятора// mass of efferent regulator orb
        public double F; //суммарная нагрузка на вал паровой машины//total load on steam-engine axle
        public double g; //ускорение свободного падения// acceleration of gravity
        public double J; //момент инерции маховика паровой машины//moment of inertia steem-engine flywheel
        public double n; //передаточное число//gear-ratio
        public double b; //коэффициент трения//constant of friction
        public double k; //константа регуляции//regulation constant

        public regulatorConditions()
        {
            m = 1;
            F = 1;
            g = 10;
            J = 1;
            n = 0.5;
            b = 0.25;
            k = 2;
            w = null;
            f = null;
            df = null;
            h = null;
            s = null;

        }
        public regulatorConditions(double m, double F, double g, double J, double n, double b, double k, double w0, double f0, double df0, double h0, int maxStep)
        {
            this.m = m;
            this.F = F;
            this.g = g;
            this.J = J;
            this.n = n;
            this.b = b;
            this.k = k;
            this.maxStep = maxStep;

            w = new double[maxStep + 2];
            w[0] = w0;
            f = new double[maxStep + 2];
            f[0] = f0;
            df = new double[maxStep + 2];
            df[0] = df0;

            h = new double[maxStep + 2];
            h[0] = h0;

            s = new double[maxStep + 2];


        }

        public double[] Eqw(double wSt, double fSt, double dfSt)  //вщзвращает массив производных w, fi, dfi
        {
            double[] res = new double[3];

            res[0] = (k * Math.Cos(fSt) - F)/J;
            res[1] = dfSt;
            res[2] = n * n * wSt * wSt * Math.Sin(fSt) * Math.Cos(fSt) - 10 * Math.Sin(fSt) - b * dfSt/m;

            return res;
        }

        public double[] RK(double h, double wSt, double fSt, double dfSt) // w, f, df
        {
            double[] res = new double[3];

            double[] k1 = new double[3];
            k1[0] = Eqw(wSt, fSt, dfSt)[0];
            k1[1] = Eqw(wSt, fSt, dfSt)[1];
            k1[2] = Eqw(wSt, fSt, dfSt)[2];

            double[] k2 = new double[3];
            k2[0] = Eqw(wSt + k1[0] * h / 2, fSt + k1[1] * h / 2, dfSt + k1[2] * h / 2)[0];
            k2[1] = Eqw(wSt + k1[0] * h / 2, fSt + k1[1] * h / 2, dfSt + k1[2] * h / 2)[1];
            k2[2] = Eqw(wSt + k1[0] * h / 2, fSt + k1[1] * h / 2, dfSt + k1[2] * h / 2)[2];

            double[] k3 = new double[3];
            k3[0] = Eqw(wSt + k2[0] * h / 2, fSt + k2[1] * h / 2, dfSt + k2[2] * h / 2)[0];
            k3[1] = Eqw(wSt + k2[0] * h / 2, fSt + k2[1] * h / 2, dfSt + k2[2] * h / 2)[1];
            k3[2] = Eqw(wSt + k2[0] * h / 2, fSt + k2[1] * h / 2, dfSt + k2[2] * h / 2)[2];

            double[] k4 = new double[3];
            k4[0] = Eqw(wSt + k3[0] * h, fSt + k3[1] * h, dfSt + k3[2] * h)[0];
            k4[1] = Eqw(wSt + k3[0] * h, fSt + k3[1] * h, dfSt + k3[2] * h)[1];
            k4[2] = Eqw(wSt + k3[0] * h, fSt + k3[1] * h, dfSt + k3[2] * h)[2];

            res[0] = wSt + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0])/6;
            res[1] = fSt + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])/6;
            res[2] = dfSt + h * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2])/6;

            return res;
        }

        public void Task(double eps)
        {
            int i = 1;
            double[] wTmp = new double[2];
            double[] fTmp = new double[2];
            double[] dfTmp = new double[2];

            while (i <= maxStep)
            {



                w[i] = RK(h[i - 1], w[i - 1], f[i - 1], df[i - 1])[0];
                f[i] = RK(h[i - 1], w[i - 1], f[i - 1], df[i - 1])[1];
                df[i] = RK(h[i - 1], w[i - 1], f[i - 1], df[i - 1])[2];

                wTmp[0] = RK(h[i - 1] / 2, w[i - 1], f[i - 1], df[i - 1])[0];
                fTmp[0] = RK(h[i - 1] / 2, w[i - 1], f[i - 1], df[i - 1])[1];
                dfTmp[0] = RK(h[i - 1] / 2, w[i - 1], f[i - 1], df[i - 1])[2];

                wTmp[1] = RK(h[i - 1] / 2, wTmp[0], fTmp[0], dfTmp[0])[0];
                fTmp[1] = RK(h[i - 1] / 2, wTmp[0], fTmp[0], dfTmp[0])[1];
                dfTmp[1] = RK(h[i - 1] / 2, wTmp[0], fTmp[0], dfTmp[0])[2];

                s[i - 1] = Math.Abs(wTmp[1] - w[i]);
                if (s[i - 1] < Math.Abs(fTmp[1] - f[i])) s[i - 1] = Math.Abs(fTmp[1] - f[i]);
                if (s[i - 1] < Math.Abs(dfTmp[1] - df[i])) s[i - 1] = Math.Abs(dfTmp[1] - df[i]);
                s[i - 1] = s[i - 1] / Math.Pow(2.0, 5);

                if (Math.Abs(s[i - 1]) < eps / Math.Pow(2.0, 5))
                {
                    h[i] = 2 * h[i - 1];
                    i++;
                }
                else
                {
                    if (Math.Abs(s[i - 1]) > eps)
                    {
                        h[i - 1] = h[i - 1] / 2;
                        
                    }
                    else
                    {
                        h[i] = h[i - 1];
                        
                        i++;
                    }
                }



            }
        }

    }
}
